Week 12 · Classification and Clustering

Set-up

We will be using the developing tidymodels framework this week for integrating with the different machine-learning libraries in a consistent manner. You can install this package from the usual RStudio Tools menu. All of the other tools that are strictly necessary for clustering are available in base R. For full flexibility, however, the ggdendro, protoclust, and heatmaply packages are recommended. If you want to explore further possibilities, look at the cluster package.

As you work through this lab, you will occasionally get error messages asking you to install other packages, too. Install whatever R asks for from the Tools menu, and then try running the chunk again.

library(tidyverse)
library(tidymodels)
library(protoclust)
library(ggdendro)
library(heatmaply)
library(spotifyr)
library(compmus)
source('spotify.R')

Classification

In order to demonstrate some of the principles of classification, we will try to identify some of the features that Spotify uses to designate playlists as ‘workout’ playlists. For a full analysis, we would need to delve deeper, but let’s start with a comparison of three playlists: Indie Pop, Indie Party, and Indie Workout. For speed, this example will work with only the first 20 songs from each playlist, but you should feel free to use more if your computer can handle it.

After you have this section of the notebook working, try using other combinations of Spotify workout playlists with similarly-named non-workout playlists.

pop <- 
    get_playlist_audio_features('spotify', '37i9dQZF1DWWEcRhUVtL8n') %>% 
    slice(1:20) %>% 
    add_audio_analysis
party <- 
    get_playlist_audio_features('spotify', '37i9dQZF1DWTujiC7wfofZ') %>% 
    slice(1:20) %>% 
    add_audio_analysis
workout <- 
    get_playlist_audio_features('spotify', '37i9dQZF1DXaRL7xbcDl7X') %>% 
    slice(1:20) %>% 
    add_audio_analysis

As you think about this lab session – and your portfolio – think about the four kinds of validity that Sturm and Wiggins discussed in our reading for last week. Do these projects have:

• Statistical validity [somewhat beyond the scope of this course]?
• Content validity?
• Internal validity?
• External validity?

We bind the three playlists together using the trick from Week 7, transpose the chroma vectors to a common tonic using the compmus_c_transpose function, and then summarise the vectors like we did when generating chromagrams and cepstrograms. Again, Aitchison’s clr transformation can help with chroma.

indie <- 
    pop %>% mutate(playlist = "Indie Pop") %>% 
    bind_rows(
        party %>% mutate(playlist = "Indie Party"),
        workout %>% mutate(playlist = "Indie Workout")) %>% 
    mutate(playlist = factor(playlist)) %>% 
    mutate(
        segments = 
            map2(segments, key, compmus_c_transpose)) %>% 
    mutate(
        pitches = 
            map(segments, 
                compmus_summarise, pitches, 
                method = 'mean', norm = 'manhattan'),
        timbre =
            map(
                segments,
                compmus_summarise, timbre,
                method = 'mean')) %>% 
    mutate(pitches = map(pitches, compmus_normalise, 'clr')) %>% 
    mutate_at(vars(pitches, timbre), map, bind_rows) %>% 
    unnest(cols = c(pitches, timbre))

Pre-processing

In the tidyverse approach, we can preprocess data with a recipe specifying what we are predicting and what variables we think might be useful for that prediction. Then we use step functions to do any data cleaning (usually centering and scaling, but step_range is a viable alternative that squeezes everything to be between 0 and 1). Finally we prep and juice the data.

indie_class <- 
    recipe(playlist ~
               danceability +
               energy +
               loudness +
               speechiness +
               acousticness +
               instrumentalness +
               liveness +
               valence +
               tempo +
               duration +
               C + `C#|Db` + D + `D#|Eb` +
               E + `F` + `F#|Gb` + G +
               `G#|Ab` + A + `A#|Bb` + B +
               c01 + c02 + c03 + c04 + c05 + c06 +
               c07 + c08 + c09 + c10 + c11 + c12,
           data = indie) %>% 
    step_center(all_predictors()) %>%
    step_scale(all_predictors()) %>%
    # step_range(all_predictors()) %>% 
    prep(indie) %>% 
    juice

Cross-Validation

The vfold_cv function sets up cross-validation. We will use 5-fold cross-validation here in the interest of speed, but 10-fold cross-validation is more typical.

indie_cv <- indie_class %>% vfold_cv(5)

Classification Algorithms

Your DataCamp tutorials this week introduced four classical algorithms for classification: \(k\)-nearest neighbour, naive Bayes, logistic regression, and decision trees. Other than naive Bayes, all of them can be implemented more simply in tidymodels. In order to use cross-validation, however, we need to write some local helper functions to fit the classifier on the training sets, predict the labels for the test/validation sets, and bind the results to the original data.

\(k\)-Nearest Neighbour

A \(k\)-nearest neighbour classifier often works just fine with only one neighbour. It is very sensitive to the choice of features, however. Let’s check the performance as a baseline and come back to it later.

indie_knn <- 
    nearest_neighbor(mode = 'classification', neighbors = 1) %>% 
    set_engine('kknn')
predict_knn <- function(split)
    fit(indie_knn, playlist ~ ., data = analysis(split)) %>% 
    predict(assessment(split), type = 'class') %>%
    bind_cols(assessment(split))

After a little awkwardness with cross-validation, we can use conf_mat to get a confusion matrix.

indie_cv %>% 
    mutate(pred = map(splits, predict_knn)) %>% unnest(pred) %>% 
    conf_mat(truth = playlist, estimate = .pred_class)

               Truth
Prediction      Indie Party Indie Pop Indie Workout
  Indie Party             8         2             5
  Indie Pop               5         7            10
  Indie Workout           7        11             5

These matrices autoplot in two forms.

indie_cv %>% 
    mutate(pred = map(splits, predict_knn)) %>% unnest(pred) %>% 
    conf_mat(truth = playlist, estimate = .pred_class) %>% 
    autoplot(type = 'mosaic')

indie_cv %>% 
    mutate(pred = map(splits, predict_knn)) %>% unnest(pred) %>% 
    conf_mat(truth = playlist, estimate = .pred_class) %>% 
    autoplot(type = 'heatmap')

We can also compute statistics like accuracy, Cohen’s kappa, or the J-measure. (Cohen’s kappa and the J-measure are two popular methods to account for how well a classifier would do just by random chance; they also range from 0 to 1.)

indie_cv %>% 
    mutate(pred = map(splits, predict_knn)) %>% unnest(pred) %>% 
    metric_set(accuracy, kap, j_index)(truth = playlist, estimate = .pred_class)

Logistic and Multinomial Regression

In the two-class case, we use logistic regression, but beware if you have more than two classes! R will just build a classifier for the first two without warning.

indie_logistic <- 
    logistic_reg(mode = 'classification') %>% 
    set_engine('glm')
predict_logistic <- function(split)
    fit(indie_logistic, playlist ~ ., data = analysis(split)) %>% 
    predict(assessment(split), type = 'class') %>%
    bind_cols(assessment(split))

With three or more classes, we need multinomial regression instead. You can adjust the penalty parameter if you are feeling adventurous.

indie_multinom <- 
    multinom_reg(mode = 'classification', penalty = 0.1) %>% 
    set_engine('glmnet')
predict_multinom <- function(split)
    fit(indie_multinom, playlist ~ ., data = analysis(split)) %>% 
    predict(assessment(split), type = 'class') %>%
    bind_cols(assessment(split))

It is not a strong classifier for this problem.

indie_cv %>% 
    mutate(pred = map(splits, predict_multinom)) %>% unnest(pred) %>% 
    metric_set(accuracy, kap, j_index)(truth = playlist, estimate = .pred_class)

We can look at the most important features in the model by using the coef method.

indie_class %>%
    fit(indie_multinom, playlist ~ ., data = .) %>%
    pluck('fit') %>%
    coef(s = 0)

$`Indie Party`
35 x 1 sparse Matrix of class "dgCMatrix"
                           1
(Intercept)      -1.78285212
danceability     -2.53638971
energy            .         
loudness          .         
speechiness      -0.60388521
acousticness     -1.13132203
instrumentalness  .         
liveness          0.88748099
valence           7.65925707
tempo             1.74806561
duration          0.90272650
C                -0.85239128
`C#|Db`           4.90362948
D                 0.53885089
`D#|Eb`           .         
E                 .         
F                 2.36474010
`F#|Gb`          -2.76928962
G                 .         
`G#|Ab`           .         
A                 0.05896787
`A#|Bb`           .         
B                -0.97831090
c01               .         
c02               .         
c03               .         
c04               .         
c05              -0.00371209
c06               .         
c07               .         
c08               .         
c09              -0.57114900
c10               4.69233868
c11              -0.44930256
c12               .         

$`Indie Pop`
35 x 1 sparse Matrix of class "dgCMatrix"
                           1
(Intercept)       3.11390092
danceability      3.45026211
energy            .         
loudness          1.24589231
speechiness       4.71954826
acousticness      .         
instrumentalness -0.42555636
liveness          .         
valence          -3.80650180
tempo             .         
duration         -3.79951231
C                 .         
`C#|Db`          -4.96323745
D                -0.42936912
`D#|Eb`           .         
E                 3.81249352
F                 .         
`F#|Gb`           .         
G                 .         
`G#|Ab`           .         
A                -2.86818513
`A#|Bb`           .         
B                 0.01899732
c01               .         
c02               .         
c03               .         
c04               .         
c05               9.42624955
c06               .         
c07              -0.57712744
c08               3.45576116
c09               .         
c10               .         
c11               .         
c12               .         

$`Indie Workout`
35 x 1 sparse Matrix of class "dgCMatrix"
                             1
(Intercept)       -1.331048804
danceability       .          
energy            -0.008703007
loudness          -0.462533428
speechiness        .          
acousticness       2.343274705
instrumentalness   1.169862384
liveness         -12.314715746
valence            .          
tempo             -1.892540056
duration           .          
C                  .          
`C#|Db`            .          
D                  .          
`D#|Eb`            2.706392941
E                  .          
F                 -0.697028182
`F#|Gb`            0.178237512
G                  .          
`G#|Ab`            .          
A                  .          
`A#|Bb`            .          
B                  .          
c01                .          
c02               -2.519830116
c03                .          
c04                5.116776480
c05                .          
c06               -0.022005255
c07                6.832908719
c08               -2.452948641
c09                .          
c10                .          
c11                0.704548663
c12                .          

Decision Trees

Decision trees are nicely intuitive, and perform somewhat better here.

indie_tree <- 
    decision_tree(mode = 'classification') %>%
    set_engine('C5.0')
predict_tree <- function(split)
    fit(indie_tree, playlist ~ ., data = analysis(split)) %>% 
    predict(assessment(split), type = 'class') %>%
    bind_cols(assessment(split))

indie_cv %>% 
    mutate(pred = map(splits, predict_tree)) %>% unnest(pred) %>% 
    metric_set(accuracy, kap, j_index)(truth = playlist, estimate = .pred_class)

We can look at the whole tree with the summary command. Be careful not to read too much into the actual numerical values, however: remember that the features were standardised before we started classification. Without cross-validation, the algorithm looks much better from the summary than it actually was in practice, but we can still see that timbre features are important and chroma features probably aren’t.

indie_class %>% 
    fit(indie_tree, playlist ~ ., data = .) %>% 
    pluck('fit') %>%
    summary

Call:
C5.0.default(x = x, y = y, trials = 1, control = C50::C5.0Control(minCases = 2, sample = 0))


C5.0 [Release 2.07 GPL Edition]     Wed Mar 18 01:24:26 2020
-------------------------------

Class specified by attribute `outcome'

Read 60 cases (35 attributes) from undefined.data

Decision tree:

acousticness > 1.061204: Indie Pop (12/1)
acousticness <= 1.061204:
:...acousticness > 0.1623779: Indie Workout (7)
    acousticness <= 0.1623779:
    :...c10 <= -0.8386642:
        :...speechiness <= -0.3888929: Indie Workout (4)
        :   speechiness > -0.3888929: Indie Pop (2)
        c10 > -0.8386642:
        :...c09 > 1.307889: Indie Pop (3)
            c09 <= 1.307889:
            :...c03 <= -2.067993: Indie Workout (2)
                c03 > -2.067993:
                :...F > 0.3436538:
                    :...energy <= -0.1647938: Indie Pop (3)
                    :   energy > -0.1647938: Indie Party (5/1)
                    F <= 0.3436538:
                    :...`A#\|Bb` > 1.649611: Indie Workout (2)
                        `A#\|Bb` <= 1.649611:
                        :...c07 > 1.513791: Indie Workout (2)
                            c07 <= 1.513791:
                            :...valence > -0.7504143: Indie Party (12)
                                valence <= -0.7504143:
                                :...c08 <= 0.1367396: Indie Workout (2)
                                    c08 > 0.1367396: Indie Party (4)


Evaluation on training data (60 cases):

        Decision Tree   
      ----------------  
      Size      Errors  

        13    2( 3.3%)   <<


       (a)   (b)   (c)    <-classified as
      ----  ----  ----
        20                (a): class Indie Party
         1    19          (b): class Indie Pop
               1    19    (c): class Indie Workout


    Attribute usage:

    100.00% acousticness
     68.33% c10
     58.33% c09
     53.33% c03
     50.00% F
     36.67% `A#\|Bb`
     33.33% c07
     30.00% valence
     13.33% energy
     10.00% speechiness
     10.00% c08


Time: 0.0 secs

Random Forests

indie_forest <- 
    rand_forest(mode = 'classification') %>% 
    set_engine('randomForest')
predict_forest <- function(split)
    fit(indie_forest, playlist ~ ., data = analysis(split)) %>% 
    predict(assessment(split), type = 'class') %>%
    bind_cols(assessment(split))

indie_cv %>% 
    mutate(pred = map(splits, predict_forest)) %>% 
    unnest(pred) %>% 
    metric_set(accuracy, kap, j_index)(truth = playlist, estimate = .pred_class)

Random forests give us the best-quality ranking of feature importance, and we can plot it with randomForest::varImpPlot. Again, it is clear that timbre, specifically Component 1 (power) and Component 11, is important. Note that because random forests are indeed random, the accuracy and feature rankings will vary (slightly) every time you re-run the code.

indie_class %>% 
    fit(indie_forest, playlist ~ ., data = .) %>% 
    pluck('fit') %>% 
    randomForest::varImpPlot()

Feature Selection

Let’s try \(k\)-NN again with just the top features. We see much better results.

predict_knn_reduced <- function(split)
    fit(
        indie_knn, 
        playlist ~ c01 + liveness + acousticness + c02 + energy, 
        data = analysis(split)) %>% 
    predict(assessment(split), type = 'class') %>%
    bind_cols(assessment(split))
indie_cv %>% 
    mutate(pred = map(splits, predict_knn_reduced)) %>% unnest(pred) %>% 
    metric_set(accuracy, kap, j_index)(truth = playlist, estimate = .pred_class)

indie_cv %>% 
    mutate(pred = map(splits, predict_knn_reduced)) %>% unnest(pred) %>% 
    conf_mat(truth = playlist, estimate = .pred_class) %>% 
    autoplot(type = 'mosaic')

Armed with this feature set, perhaps we can make a better plot. It’s clear that the workout list has fewer live tracks, and that the party playlist is somewhat louder and higher on Components 1 and 2 than the pop list.

indie %>%
    ggplot(aes(x = c01, y = c02, colour = playlist, size = liveness)) +
    geom_point(alpha = 0.8) +
    scale_color_brewer(type = 'qual', palette = 'Accent') +
    labs(
        x = 'Timbre Component 1', 
        y = 'Timbre Component 2', 
        size = 'Liveness', 
        colour = 'Playlist'
    )

Can you get better performance by using more or fewer of the top features from the random forest?

Clustering (Optional)

The Bibliothèque nationale de France (BnF) makes a large portion of its music collection available on Spotify, including an eclectic collection of curated playlists. The defining musical characteristics of these playlists are sometimes unclear: for example, they have a Halloween playlist. Perhaps clustering can help us organise and describe what kinds of musical selections make it into the BnF’s playlist.

We begin by loading the playlist and summarising the pitch and timbre features, just like last week. Note that, also like last week, we use compmus_c_transpose to transpose the chroma features so that – depending on the accuracy of Spotify’s key estimation – we can interpret them as if every piece were in C major or C minor. Although this example includes no delta features, try adding them yourself if you are feeling comfortable with R!

halloween <- 
    get_playlist_audio_features('bnfcollection', '1vsoLSK3ArkpaIHmUaF02C') %>% 
    add_audio_analysis %>% 
    mutate(
        segments = 
            map2(segments, key, compmus_c_transpose)) %>% 
    mutate(
        pitches = 
            map(segments, 
                compmus_summarise, pitches, 
                method = 'mean', norm = 'manhattan'),
        timbre =
            map(
                segments,
                compmus_summarise, timbre,
                method = 'mean')) %>% 
    mutate(pitches = map(pitches, compmus_normalise, 'clr')) %>% 
    mutate_at(vars(pitches, timbre), map, bind_rows) %>% 
    unnest(cols = c(pitches, timbre))

Pre-processing

Remember that in the tidyverse approach, we can preprocess data with a recipe. In this case, instead of a label that we want to predict, we start with a label that will make the cluster plots readable. For most projects, the track name will be the best choice (although feel free to experiment with others). The code below uses str_trunc to clip the track name to a maximum of 20 characters, again in order to improve readability. The other change from last week is column_to_rownames, which is necessary for the plot labels to appear correctly.

Last week we also discussed that although standardising variables with step_center to make the mean 0 and step_scale to make the standard deviation 1 is the most common approach, sometimes step_range is a better alternative, which squashes or stretches every features so that it ranges from 0 to 1. For most classification algorithms, the difference is small; for clustering, the differences can be more noticable. It’s wise to try both.

halloween_juice <- 
    recipe(track.name ~
               danceability +
               energy +
               loudness +
               speechiness +
               acousticness +
               instrumentalness +
               liveness +
               valence +
               tempo +
               duration +
               C + `C#|Db` + D + `D#|Eb` +
               E + `F` + `F#|Gb` + G +
               `G#|Ab` + A + `A#|Bb` + B +
               c01 + c02 + c03 + c04 + c05 + c06 +
               c07 + c08 + c09 + c10 + c11 + c12,
           data = halloween) %>% 
    step_center(all_predictors()) %>%
    step_scale(all_predictors()) %>%
    # step_range(all_predictors()) %>% 
    prep(halloween %>% mutate(track.name = str_trunc(track.name, 20))) %>% 
    juice %>% 
    column_to_rownames('track.name')

Computing distances

When using step_center and step_scale, then the Euclidean distance is usual. When using step_range, then the Manhattan distance is also a good choice: this combination is known as Gower’s distance and has a long history in clustering.

After you have this section of the notebook working with Euclidean distance, try modifying it to use Gower’s distance.

halloween_dist <- dist(halloween_juice, method = 'euclidean')

Hierarchical clustering

As you learned in your DataCamp exercises this week, there are three primary types of linkage: single, average, and complete. Usually average or complete give the best results. We can use the ggendrogram function to make a more standardised plot of the results.

hclust(halloween_dist, method = 'single') %>% dendro_data %>% ggdendrogram

A more recent – and often superior – linkage function is minimax linkage, available in the protoclust package. It is more akin to \(k\)-means: at each step, it chooses an ideal centroid for every cluster such that the maximum distance between centroids and all members of their respective clusters is as small as possible.

protoclust(halloween_dist) %>% dendro_data %>% ggdendrogram

Try all four of these linkages. Which one looks the best? Which one sounds the best (when you listen to the tracks on Spotify)? Can you guess which features are separating the clusters?

k-Means

Unlike hierarchical clustering, k-means clustering returns a different results every time. Nonetheless, it can be a useful reality check on the stability of the clusters from hierarchical clustering.

Try different numbers of clusters and see which results are the most stable.

kmeans(halloween_juice, 4)

K-means clustering with 4 clusters of sizes 5, 2, 7, 6

Cluster means:
  danceability      energy    loudness speechiness acousticness instrumentalness    liveness
1   -0.9168866 -0.09565648  0.04274727 -0.62837629   -0.2771348        0.4410513  0.56188821
2    0.5081080  0.06821405  0.19540272 -0.08024282   -0.3408504        0.4588850 -0.89728005
3    0.8866356  0.93909123  0.84609924  0.85021922   -0.3662239       -0.2242270 -0.01779592
4   -0.4397054 -1.03863072 -1.08787274 -0.44152791    0.7718237       -0.2589062 -0.14838492
     valence       tempo   duration            C      C#|Db          D       D#|Eb          E
1 -0.4536494 -0.36865124 -0.2179715  0.258484457 -0.6833355  0.4248468 -0.21559098  1.1969673
2  0.5451643  0.19468869 -0.5870669 -0.875041414  1.4179889  1.5859285 -1.44260505 -1.3583614
3  0.8768736  0.17842765 -0.2982245  0.062311297 -0.1287512 -0.3072340  0.04728937 -0.3249789
4 -0.8266995  0.03414755  0.7252605  0.003580244  0.2469930 -0.5242421  0.60535657 -0.1655436
            F      F#|Gb           G       G#|Ab          A      A#|Bb          B        c01
1 -0.18647403 -0.6425998  0.08079254 -0.94006032  0.4212267 -0.9058901  0.8885227  0.7812260
2 -0.67457284  2.2071822  1.30025415  0.05420713 -0.4323913 -0.8310957 -0.5304488 -0.3034831
3  0.06016135  0.1576683 -0.50878092  0.13519034  0.2728931  0.4500771  0.1026522  0.4571862
4  0.31006440 -0.3841739  0.09283258  0.60759250 -0.5252671  0.5068503 -0.6833802 -1.0832445
         c02         c03         c04        c05        c06         c07        c08        c09
1  0.9616765  0.99474968 -0.04951795  0.2721605 -0.5969236 -0.54702365 -0.5949745  0.1906642
2 -0.1844395 -1.32821483  0.15230482 -0.6988863 -0.1446833  0.29754667 -0.6677161 -0.3762457
3  0.3033856 -0.40324528  0.24925101 -0.1033827  0.5991044  0.35084069  0.3799168 -0.1968883
4 -1.0938671  0.08423303 -0.30029616  0.1267749 -0.1532911 -0.05264332  0.2751478  0.1962314
         c10        c11        c12
1 -0.2892562  0.8331551  0.5543391
2  0.6902069  0.1415686  0.3740885
3  0.4170187 -0.7718281 -0.3822938
4 -0.4755440  0.1589806 -0.1406360

Clustering vector:
I Put a Spell on You      Close Your Eyes           Evil Woman         Time to Kill 
                   3                    4                    3                    2 
      Evil Bad Woman The Princess of Evil      'Round Midnight Tana's Theme - Fr... 
                   4                    4                    4                    1 
          Evil Blues  Someone Is Watching Violin Sonata in ...         Cannibal Pot 
                   3                    1                    1                    3 
           Red Devil You're Not Living...         Little Demon Devil at 4 O'Cloc... 
                   2                    1                    3                    1 
      Old Devil Moon  Up Jumped the Devil Beetween the Devi...          Devil Woman 
                   3                    3                    4                    4 

Within cluster sum of squares by cluster:
[1] 126.46108  38.93083 137.62738 126.55590
 (between_SS / total_SS =  33.5 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"   
[7] "size"         "iter"         "ifault"      

Heatmaps

Especially for storyboards, it can be helpful to visualise hierarchical clusterings along with heatmaps of feature values. We can do that with heatmaply. Although the interactive heatmaps are flashly, think carefully when deciding whether this representation is more helpful for your storyboard than the simpler dendrograms above.

grDevices::dev.size("px")

[1] 525.0000 324.4747

heatmaply(
    halloween_juice,
    hclustfun = hclust,
    # hclustfun = protoclust,
    # Comment out the hclust_method line when using protoclust.
    hclust_method = 'average',
    dist_method = 'euclidean')

Which features seem to be the most and least useful for the clustering? What happens if you re-run this section of the notebook using only the best features?